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ABSTRACT 

The quasi-steady structure of super-critical accretion flows around a black hole is studied based 
on the two-dimensional radiation-hydrodynamical (2D-RHD) simulations. The super-critical flow is 
composed of two parts: the disk region and the outflow regions above and below the disk. Within the 
disk region the circular motion as well as the patchy density structure are observed, which is caused by 
Kelvin-Helmholtz instability and probably by convection. The mass-accretion rate decreases inward, 
roughly in proportion to the radius, and the remaining part of the disk material leaves the disk to 
form outflow because of strong radiation pressure force. We confirm that photon trapping plays an 
important role within the disk. Thus, matter can fall onto the black hole at a rate exceeding the 
Eddington rate. The emission is highly anisotropic and moderately collimated so that the apparent 
luminosity can exceed the Eddington luminosity by a factor of a few in the face-on view. The mass- 
accretion rate onto the black hole increases with increase of the absorption opacity (metalicity) of the 
accreting matter. This implies that the black hole tends to grow up faster in the metal rich regions 
as in starburst galaxies or star-forming regions. 

Subject headings: accretion: accretion disks — black hole physics — hydrodynamics — radiative 
transfer 



1. INTRODUCTION 

It is widely believed that accretion flow onto black 
holes drives major activities of astrophysical black holes, 
such as active galactic nuclei (AGNs), Galactic black 
hole candidates (BHCs), and possibly gamma-ray bursts 
(GRBs). It is also a common belief that the basic ac- 
cretion processes and radiation properties can well be 
described by the standard-disk model by Shakura & Sun- 
yaev (1973). On the other hand, the observational facts 
which cannot fit the standard-disk picture are also being 
accumulated recently. Good examples are intense high- 
energy emission from black holes, which indicates the 
presence of very hot plasmas around black holes with 
temperature of T ~ 10 9 K. This leads to the ideas of 
accretion disk corona and/or radiatively inefficient flow 
(RIAF). 

The standard-disk picture breaks down not only in the 
low-luminosity regimes but also in the high-luminosity 
regimes, in which mass-accretion rates, M , becomes com- 
parable to or exceeds the critical mass-accretion rate, 
M CT it = Lb/c 2 , with Le being the Eddington lumi- 
nosity. Apparently, they look similar to those of the 
standard-type disks, since the disk remains to be opti- 
cally thick and thus emits blackbody-like emission, just 
as the standard-type disks do. Flow dynamics is, how- 
ever, distinct. 

What makes the super-critical accretion flows distinct 
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from the standard-disk-type flow is the presence of pho- 
ton trapping (Bcgelman 1978). Photon trapping occurs, 
when photon diffusion time, time for photons to travel 
from the equatorial plane to the surface, exceeds the ac- 
cretion timescale. Under such circumstances photons 
generated via the viscous process are advected inward 
with gas flow without being able to go out from the sur- 
face immediately. We can define the trapping radius, 
inside which photon trapping is significant; 
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where to is the mass-accretion rate normalized by the 
critical mass-accretion rate, H is the half thickness of 
the disk, r is the radius, and r g {= 2GM/c 2 ) is the 
Schwarzschild radius with M being the black hole mass 
(Ohsuga et al. 2002, hereafter referred to as Paper I, see 
also Begelman 1978 for the spherical case). 

It is claimed that photon trapping effects were incor- 
porated in the slim-disk model by Abramowicz et al. 
(1988; for a review, see Kato, Fukue, & Mineshige 1998, 
section 10.1). We have shown previously, however, that 
the slim-disk model does not fully consider the photon 
trapping effects, since it is a radially one-dimensional 
model, whereas the photon trapping is basically a multi- 
dimensional effect (see Paper I). We thus need at least 
two-dimensional treatment. Indeed, we have demon- 
strated in Paper I that the photon trapping is grossly 
underestimated in the slim-disk model. 

Let us recall what lacks in the slim-disk model more 
explicitly. The equation of energy balance, Q V j S = 
Qrad + Qadv, is solved in the slim-disk model, where Q v i s , 
Q ra d, and Q a dv are the viscous heating, radiative cooling, 
and advective cooling rates per unit surface, respectively. 
The problem resides in that the radiative cooling is eval- 
uated under the usual diffusion approximation (in the 
vertical direction). This approximation may be justified, 
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if radial inflow of gas is totally negligible so that photons 
can mainly diffuse in the vertical direction. However, this 
is not the case, when the diffusion timescale is longer than 
the accretion timescale. This leads to over-estimation of 
(5 ra d and, hence, under-estimations of Qadv, compared 
with the correct value. More importantly, photon trap- 
ping modifies spectral energy distribution (SED), as well. 
We have found that large photon trapping yields spectral 
softening, because hard photons which are created deep 
inside the disk are more effectively trapped than soft pho- 
tons (Ohsuga, Mineshige, & Watarai 2003). Since our 
previous formulations are only partly multi-dimensional, 
we, as a next step, need to perform fully two-dimensional 
analysis of the super-critical accretion flows. This is a 
major motivation of the present study. 

When considering photon trapping effects, we should 
also pay attention to the fact that the super-critical ac- 
cretion flow becomes geometrically thick. The multi- 
dimensional gas motion, such as convective or large-scale 
circulation, might occur. Further, strong outflow might 
also be generated at the disk surface via radiation pres- 
sure force. Such complex flow motion will influence ra- 
diation energy distribution through the advective energy 
transport, which, in turn, affects the flow motion via ra- 
diation pressure force. We need to carefully solve such 
strong coupling between radiation and matter. 

Two-dimensional radiation-hydrodynamical (2D- 
RHD) simulations of accretion disks were initiated by 
Eggum, Coroniti, & Katz (1987), who assumed the 
equilibrium between gas and radiation. The improved 
simulations, in which the energy of gas and radiation 
are separately treated, were performed by Kley (1989), 
Okuda, Fujita, & Sakashita (1997), Fujita & Okuda 
(1998), Kley & Lin (1999), and Okuda & Fujita (2000). 

The simulations of super-critical flows around black 
holes were pioneered by Eggum, Coroniti, & Katz (1988), 
which again assumed the equilibrium between gas and 
radiation, and was improved by Okuda (2002). Eggum 
et al. (1988) showed that mass accretion onto the black 
hole occurs at a super-critical rate, although the mass- 
accretion as well as mass-outflow rate were still variable 
in their simulations; that is, quasi-steady state had not 
been achieved by their simulations. This is the same 
also in the simulations by Okuda (2002), in which the 
resultant luminosity slowly decreases with time. He also 
found that the mass-accretion rate is sub-critical around 
the black hole, although the mass is injected from the 
outer boundary at a super-critical rate. Note that the 
calculation times of these simulations were only 0.6 s 
and 1.6 s (in physical unit), respectively, for the black 
hole mass of 10M Q . Notably, these are shorter than the 
viscous timescale, 
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with a being the viscosity parameter. More recently, 
Okuda et al. (2005) performed long-term 2D-RHD calcu- 
lations of the super-critical accretion. The luminosity as 
well as the mass-accretion rate seems to be quasi-steady 
in their simulations, although the flow structure is not 
steady yet. Since the sum of the mass-accretion rate at 



the inner boundary and mass-outflow rate at the outer 
boundary is much smaller than the mass-input rate at 
the disk boundary, the mass within the computational 
domain continues to increase. In addition, the photon- 
trapping did not appear in their simulations, whereas the 
mass-accretion rate exceeds the critical value. 

To summarize, despite the interesting simulations be- 
ing made so far by two groups the quasi-steady structure 
of the super-critical accretion flows still remains to be an 
open issue. Further, interesting issues related to super- 
critical flow, including the photon-trapping effects, the 
dependence of the observed luminosity on various view- 
ing angles, and the effects of metalicity on the flow struc- 
ture, have not been investigated previously. This is a 
motivation of the present study. 

Here, we report for the first time the quasi-steady 
structure of the super-critical disk accretion flows around 
black holes, which were revealed by 2D-RHD simulations. 
Through the present simulations we mainly aim at under- 
standing the dynamics of the viscous flow in the vicinity 
of the black hole of r < 100r g . Basic equations and our 

model are explained in §2, and the numerical methods 
are described in §3. We will then display the quasi-steady 
flow structure and study the photon trapping effects in 
the simulated flow in §4. Finally, §5 and §6 are devoted 
to discussion and conclusions. 

2. BASIC EQUATIONS AND ASSUMPTIONS 

We solve the full set of RHD equations including 
the viscosity term. We use spherical polar coordinates 
(r,9,ip), where r is the radial distance, 9 is the polar 
angle, and tp is the azimuthal angle. In the present 
study, we assume the flow to be non self-gravitating, re- 
flection symmetric relative to the equatorial plane (with 
9 = 7r/2), and axisymmetric with respect to the rotation 
axis (i.e., 8/ dip — 0). We describe the gravitational field 
of the black hole in terms of pseudo-Newtonian hydrody- 
namics, in which the gravitational potential is given by 
*(r) = —GMj{r — r s ) as was introduced by Paczynski 
& Wiita (1980). The basic equations are the continuity 
equation, 

^+V-(pv) = 0, (3) 
the equations of motion, 
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+ V • (pr sin 9 ■ v v v) = r sin 9 ■ q v , (6) 



the energy equation of the gas, 
de 

— + V • (ev) = -f>V • v - AttkB + ckE + $ vis , (7) 

and the energy equation of the radiation, 
dE 



dt 



■+V-(E v) = -V-Fo-Vv : P +AttkB-ckE . (8) 
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Here, p is the gas mass density, v — (v r ,vg,v v ) is the 
velocity, p is the gas pressure, e is the internal energy 
density of the gas, B is the blackbody intensity, Eg is the 
radiation energy density, Fo is the radiation flux, Po is 
the radiation pressure tensor, k is the absorption opacity, 
q = (q r ,qg,q v ) is the viscous force, and $ v is is the viscous 
dissipative function, respectively. The radiation force, 
/rad = (fr,fe), is given by 



/rad - ~F 0, 



(9) 



where x(= K+paT/m p ) is the total opacity with <tt being 
the Thomson scattering cross-section and m p being the 
proton mass. 
As the equation of state, we use 

p=( 7 -l)e, (10) 

where 7 is the specific heat ratio. The temperature of 
the gas, T, can then be calculated from 

pkT 
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where is the Boltzmann constant and p is the mean 
molecular weight, respectively. 

To complete the set of equations we apply flux limited 
diffusion (FLD) approximation developed by Levermorc 
and Pomraning (1981). In this framework, the radiation 
flux is written as 
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with the flux limiter, A, and the radiation pressure tensor 
is expressed in terms of the radiation energy density via 



Po = fEo, 



(13) 



where f is the Eddington tensor. Here, the flux limiter 
is given by 

A = 2 + U 2 , (14) 

with using the dimensionless quantity, 1Z = 
|VI?o| / ix^o)- The components of the Eddington 
tensor are 

f=^(l-/)I+^(3/-l)rm, (15) 

where / is the Eddington factor, 

f = X + \ 2 TZ 2 , (16) 

and n is the unit vector in the direction of the radiation 
energy density gradient, 

V£ 
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This approximation holds both in the optically thick 
and thin regimes. In the optically thick limit, we find 
A — ► 1/3 and / — ► 1/3 because of 1Z — > 0. In the opti- 
cally thin limit of 1Z — > 00, on the other hand, we have 
\Fq\ = cEq. These give correct relations in the optically 
thick diffusion limit and optically thin streaming limit, 
respectively. 

For the absorption opacity, we consider the free-free 
absorption, Kg, and bound- free absorption, Kbf, 



ks + km, 
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where Kg and Kbf are given by 
Kg 

(Rybicki & Lightman 1979), and 



1.7 x 1(T 25 T- 7 / 2 [ JL 



cm 



« bf = 4.8xl0- 24 T- 7 / 2 f^) (^)cm- 1 , 



(19) 



(20) 



with Z being the metalicity (Hayashi, Hoshi, & Sugimito 
1962), respectively. 

Here, we assume that only the r</?-component of the 
viscous stress tensor, which plays important roles for the 
transport of the angular momentum and heating of the 
disk plasma, is non zero, 
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in the present study, where 77 is the dynamical viscosity 
coefficient. Then, the radial and polar components of the 
viscous force are null (q r = qe = 0), and the right hand 
side of equation JBJ is described as 
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The viscous dissipative function is given by 
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Finally, we prescribe the dynamical viscosity coefficient 
as a function of the pressure 



r\ = a- 



p + XEp 



(24) 



where fix is the Keplerian angular speed. It is modi- 
fied a prescription of the viscosity, which is proposed by 
Shakura & Sunyaev (1973). In this form, the dynamical 
viscosity coefficient is proportional to the total pressure 
because of A — * 1/3 in the optically thick regime. In the 
optically thin regime, by contrast, we find 77 = ap/f^x (or 
kinematic viscosity is v = r]/p = ac 2 /£lK with c s = \pp~f~p 
being isothermal sound velocity) , since A vanishes in this 
limit. 

Here, we need to remark that the realistic formal- 
ism about the viscosity should be investigated from the 
magneto-hydrodynamical point of view, since the domi- 
nant sources of viscosity would be chaotic magnetic fields 
and turbulence in the gas flow (e.g., Machida, Hayashi, 
& Matsumoto 2000; Stone & Pringle 2001). 

3. NUMERICAL METHODS 

3.1. The code 

We numerically solve the set of RHD equations shown 
in the previous section by using an explicit-implicit finite 
difference scheme on the Eulerian grids. Our methods 
and boundary conditions are similar to those of Okuda 
(2002), but we adopt a different initial condition and 
have carried out long-term calculations in order to exam- 
ine a quasi-steady structure of the supercritical accretion 
flows. Since we assume the axisymmetry as well as the 
reflection symmetry, the computational domain can be 
restricted to one quadrant of the meridional plane. The 
domain consists of spherical shells of r- m < r < r ont and 
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< < 7r/2, where r- ln and r out are the radial coordi- 
nates of the inner and outer boundaries of the computa- 
tional domain, respectively, and is divided into N r x Ng 
grid cells, amounting to 96 x 96 meshes. The N r grid 
points in radial direction are equally spaced logarithmi- 
cally, while the Ng grids are equally distributed in such 
a way to achieve A cos (9 = 1/Ng. All the physical quan- 
tities are defined at each cell center. 

We divide the numerical procedure for the finite- 
difference equations into the following steps: [1] hydro- 
dynamical terms for ideal fluid, [2] advection term in 
the energy equation of the radiation, [3] radiation flux 
term in the radiation energy equation, [4] gas-radiation 
interaction terms in the energy equations of the radiation 
and gas, and [5] viscous terms in the momentum equa- 
tion and energy equations of the gas. Steps [1] and [2] 
are solved with the explicit method while the steps [3], 
[4] , and [5] are treated based on the implicit method. In 
the first step, we use the computational hydrodynamical 
code, the Virginia Hydrodynamics One. It is based on 
the Piecewise Parabolic Method (Colella & Woodward 
1984). The equations ©-0 except the viscous terms 
and the gas-radiation interaction terms of equation J2J 
are solved in this step. In the second step, an integral for- 
mulation is used to generate a conservative differencing 
scheme for the advection term of the equation ijHJ. The 
energy transport by the radiation flux term is solved in 
the third step. The radiation energy density is updated 
again with using the BiCGSTAB method for a matrix 
inversion. In the fourth step, we solve the gas-radiation 
interaction terms in equation Q. All the terms on the 
right hand side of equation (JSJ) except for the radiative 
flux term is also treated in this step. The radiation en- 
ergy and gas energy is advanced simultaneously. The 
method used in this step is basically the same as that 
described by Turner & Stone (2001). In the final step, 
we updates the azimuthal component of the velocity by 
solving the viscosity term in the equation JSJ with the 
Gauss- Jordan elimination for a matrix inversion. The 
viscous dissipative function is also calculated in this step. 

Throughout the present study, we assume M = 10M©, 
a = 0.1, 7 = 5/3, and p, = 0.5. The size of the compu- 
tational domain is set to be (n n ,r ou t) = (3r g , 500r g ). 
[Here, it is noted that our conclusion does not change so 
much, even if we employ n n = 2r g .] 

3.2. Boundary and Initial Conditions 

We employ the absorbing inner boundary condition so 
that the density, gas pressure, and velocity are smoothly 
dumped (i.e. Kato, Mineshige, & Shibata 2004). Here, 
it is noted that our simulation results are not sensitive to 
the inner boundary condition nor to the location of the 
inner boundary. The results does not change even if we 
use the free boundary conditions, where all the matter 
and waves can transmit freely. The radiation flux is set 
to be F£ — —cEo at the inner boundary; that is, we 
apply the condition of the optically thin limit at the inner 
boundary. 

The outer boundary at r = r out is divided into two 
parts: the disk part (8 > 0.457r) and the part above 
the disk (9 < 0.457r). Through the outer-disk bound- 
ary we assume that matter is continuously injected into 
the computational domain. We set the injected mass- 
accretion rate (mass-input rate) so as to be constant 



at the boundary. In the present study, we explore the 
cases of mi npu t = 300, 1000, and 3000, where rhi npu t 
is the mass-input rate normalized by the critical mass- 
accretion rate. The injected matter is supposed to rotate 
with sub-Keplerian speed, and it has a specific angular 
momentum corresponding to the Keplerian angular mo- 
mentum at r = 100r g , since our main purpose of the 
present study is to investigate the viscous accretion flows 
within r < 100r g . By setting the Keplerian radius (100r g ) 

to be much smaller than the radius of the outer bound- 
ary, we can prevent that the outer boundary conditions 
directly affect the accretion flow within 100r g . With such 
a large outer boundary we can reproduce complex mo- 
tions like circulation around 100r g in the disk; that is, 
the matter can transiently go out across the radius of 
100r g and return. At the outer boundary region above 
the accretion disk we use free boundary conditions and 
allow for matter to go out but not to come in. If the 
radial velocity is negative at the outermost grid, it is au- 
tomatically set to be zero. We also employ the radiation 
flux in the optically thin limit, Fq ~ cE$, at the outer 
boundary. 

With respect to the the rotation axis we assume p, p, 
v r , and Eq to be symmetric, while vg and v v are antisym- 
metric. On the equatorial plane, on the other hand, p, p, 
v r , v v , and Eq are symmetric and vg is antisymmetric. 

We start the calculations with a hot, rarefied, and 
optically-thin atmosphere. There is no cold dense disk 
in the computational domain, initially. The initial at- 
mosphere is constructed to approximately achieve hy- 
drostatic equilibrium in the radial direction; namely, its 
density profile is given by 



P = Pout exp 



fim p GM /rout 
kBTr out V r 
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(25) 



where p ut is the density at the outer boundary. We em- 
ploy p ou t = 10~ 17 gcm~ 3 and T — 10 11 K in the present 
study. Since this atmospheric gas is finally ejected out 
of the computational domain, it does not affect the re- 
sultant quasi-steady structure. 

3.3. Time Step 

The time step is restricted by the Courant-Friedrichs- 
Levi condition. We set the time step as 



At = £ min 
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where £ is a parameter, and Ar and A9 are the cell sizes 
in the radial and polar directions, respectively. Although 
we only show the results for the cases with £ — 0.4, we 
also simulated the cases with £ — 0.1 or 0.05, confirming 
that our conclusions do not alter by this change. 

4. RESULTS 

4.1. Quasi-steady Structure 

We first represent time evolution of 2D-RHD simu- 
lations. Overall evolution is divided into two distinct 
phases: the accumulation phase and the quasi-steady 
phase. 

The mass is continuously injected through the outer- 
disk boundary, and creates continuous gas inflow because 
of the gravity force by the central black hole. Since an- 
gular momentum of the injected mass is set to be equal 



around Black Holes 



5 




2D-RHD simulations. The critical time (~ 7s) separat- 
ing these two phases roughly coincides with the viscous 
timescale [equation ©]. It is natural that the quasi- 
steady structure is obtained within 10 sec, which is the 
viscous timescale at the Keplerian radius (r ~ 100r g ), 
although the size of the computational domain (500r g ) 
is much larger. This is because that the injected mat- 
ter accretes from the outer boundary to the Keplerian 
radius with the free-fall velocity so that the evolution- 
ary timescale of the outer regions is given by the free- 
fall timescale at r = 500r g , 1.6s (r/500r g ) 3/2 (M/10M Q ). 
This is certainly shorter than the viscous timescale at 
r = IOOtv. 



50 100 150 200 

Fig. 1. — The 2D density distribution at the elapsed time of 
t =5.9 s for M = IOMq. The adopted parameters are minput 

to the Keplerian angular momentum at r = 100r g , it 
is natural that the gas tends to accumulate around the 
regions with the radius of 100r g by degrees (see Figure 
This is the accumulation phase. Since Eggum et al. 
(1988) and Okuda (2002) started calculations with a cold 
dense disk, this phase did not appear in their simulations. 
Eventually, the viscosity starts to work so that the an- 
gular momentum of the gas can be transported outward, 
which drives inflow gas motion in a quasi-steady fashion. 
This is the quasi-steady phase. 

Such a two-step evolution is clear in Figure In 
this figure, we show the time evolution of normalized 
mass-accretion rate onto the central black hole, m = 
M/(Le/c 2 ) (thick solid curve), the luminosity, L/Le 
(thin solid curve), the viscous heating rate, L v [ s /Le (dot- 
ted curve), and the total mass contained within the com- 
putational domain, m to tai = M tota i/10 19 gcm- 3 (dashed 
curve). Here, we set the mass- input rate of mi npu t = 1000 
and the metalicity of Z = Z Q . Both of the accretion rate 
and luminosity steadily increases until t<7 s. We also 
notice that the mass- accretion rate is much smaller than 
the mass-input rate (rh.input = 1000), and that the total 
mass within the computational domain rapidly increases 
with time, indicating mass accumulation really occurring 
in this phase. Then starts the quasi-steady accretion 
phase (at t>7s), when the mass-accretion rate exceeds 
critical value, rh > 1, and all the physical quantities stay 
nearly constant. [The constant 77i to tai implies that the 
sum of the mass-accretion rate at the inner boundary 
and mass-outflow rate at the outer boundary is equal 
to the mass-input rate.] We can thus conclude that we 
could for the first time succeed in reproducing the quasi- 
steady state of the super-critical accretion flows with 
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Fig. 2. — Time evolutions of the total mass within the compu- 
tational domain (dashed curve), the mass-accretion rate onto the 
black hole normalized by Le/c 2 (thick solid curve), the luminosity 
normalized by (thin solid curve), and the total viscous heating 
rate normalized by Lj; (dotted curve), respectively, for the case 
with rhi nput = 1000 and Z = Zq. 

Let us, next, examine the quasi-steady structure in 
some details. Figure displays the cross-sectional view 
of the density distributions (with colors), overlaid with 
the velocity vectors (with arrows) at t = 10 s. Here, x- 
and y-axis are R — rsinO and z = rcosO, respectively. 
We understand with this figure that the flow structure is 
roughly divided into two regions: the disk region around 
the equatorial plane (characterized by orange color) and 
the outflow region above the inflow region (characterized 
by blue color). Roughly, the boundary is at z/R ~ 0.8; 
that is the disk is geometrically and optically thick, as 
was predicted by the slim-disk model. However, the den- 
sity distribution definitely deviates from that of the slim- 
disk model, since it is neither smooth nor plane parallel 
in the vertical direction. We can even see a number of 
cavities in this figure. The flow pattern is also complex, 
though the slim-disk model predicts the simple conver- 
gence flow. We found the prominent circular motion 
within the disk and the strong outflow which is gener- 
ated at the disk surface. The patchy structure around 
the boundary between the dense inflow region and the 
rarefied outflow region seems to be caused by the Kelvin- 
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Fig. 3. — The 2D density distribution overlaid with the velocity vectors at the elapsed time of t =10 s for M = IOMq. The adopted 



parameters are m\ 



nput 



1000 and Z - 



Helmholtz (K-H) instability (discussed later). The less 
dense gas in the outflow region penetrates into the disk 
body because of the K-H instability, thus forming the 
cavities. 

To understand the dynamics of our entire calcula- 
tion domain, we show the 2D density distribution in 
the whole region in Figure 0] We see from this figure 
that the density is larger around the equatorial plane 
(|z|<lGT g ) than at the large altitudes (|z|>10r g ) in 
the outer region, i?>100r g . This reflects that the in- 
jected mass accretes along the equatorial plane. On 
the other hand, the viscous accretion disk forms inside, 
r<100r g . We focus the viscous accretion flows in the 

present study. The situation is apparently similar to 
that studied by Chakrabarti (1996), although they were 
concerned with subcritical flow and their study was not 



multi-dimensional nor radiation-hydrodynamical simula- 
tions. 

Figure [5] show the 2D distributions of radiation energy 
density (upper left), the ratio of the radiation energy to 
the internal energy of gas (upper right), gas temperature 
(lower left), and the radial velocity normalized by the 
escape velocity (lower right), respectively, on the R — z 
plane. As shown in the upper left panel, the radiation en- 
ergy distribution roughly coincides with the gas density 
distribution. That is, the radiation energy tends to be 
larger around the equatorial plane than that around the 
rotation axis. Since the radiation energy distribution is 
smoothed due to the radiative diffusion within the disk, 
there is no cavity found in the radiation energy distribu- 
tion, which makes a marked difference from the density 
distribution (see Figure QJ. 

Radiation energy greatly exceeds the gas energy in the 
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Fig. 4. — The 2D density distribution in the whole computational domain at the elapsed time of t =10 s for M = IOMq. The adopted 
parameters are mi nput = 1000 and Z = Zq. 



entire region, including the outflow region in our simu- 
lations (see the upper right panel). This confirms that 
most of the regions is the radiation-pressure dominated 
and strong radiation pressure supports the geometrically 
thick disk and drives the outflow. 

The gas temperature distribution shown in the lower 
left panel shows relatively low temperatures in the disk 
region, compared with those in the outflow region, al- 
though the viscous heating rate is much larger in the 
former than in the latter. This can be understood, be- 
cause radiative cooling is more effective in dense region 
as a consequence of its strong density dependence. The 
gas in the disk region is heated up by the viscous heating, 
and the processed energy is effectively converted to the 
radiation energy. Therefore, gas temperature does not 
rise so much within the disk. Conversely, the gas cannot 
emit effectively in the outflow region, since both free-free 



emissivity and bound-free emissivity are more sensitive 
to the density than the gas temperature, oc p 2 T 1 / 2 . It 
can be understood by the comparison between the lower 
left and the upper right panels. The radiation tempera- 
ture (oc E^ 4 ) in the outflow region is lower than that in 
the disk region, on the contrary to the gas temperature 
profile. 

The lower right panel indicates the radial velocity nor- 
malized by the escape velocity. The gas moves toward 
the black hole or flows outward slowly in the disk re- 
gions (see the blue area). The white color indicates that 
the velocity exceeds the escape velocity in this area. It is 
found that the gas is accelerated through the radiation 
pressure and is blown away to a large distance (see also 
the upper right panel). Such a flow component will be 
identified as a strong disk wind. Note that the outflow 





Fig. 5. — The 2D distribution of the radiation energy density (upper left), the ratio of the radiation energy to the inertial energy of the 
gas (upper right), the gas temperature (lower left), and the radial velocities normalized by the escape velocity (lower right), respectively, 
at t = 10 s. The adopted parameters are rh; nput = 1000 and Z = Zq . 



presented here is distinct in nature from that known as 
a bounce jet (Chen et al. 1997), which arises because of 
a bounce of free-falling low angular momentum material, 
when it goes through the centrifugal barrier at small r. 

The outflow will also produce large absorption in the 
emergent spectra. We also notice strong velocity shear 
at the boundary between the disk region and the outflow 
region. This complex density profile around the disk sur- 
face as shown in Figure |21 is explained as a consequence 
of the K-H instability. 



The growth timescale of the K-H instability is roughly 
given by 



£kh 



1 

kv r 



Pdisk 
Pout 



1/2 



(27) 



where k is the wavenumber, while pdisk/ Pout is the den- 
sity ratio of the disk region to the outflow region at the 
disk boundary. Here, we assume the incompressible fluid 
as well as pdisk 3> Pout and neglect the gravity. Also, 
the viscosity is not considered, since r#-coniponent of 
the viscous stress tensor is set to be zero in the present 
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simulations. By setting k — 27r/10r g , v r = 0.1c, and 
Pdisk/Pout = 10, we find £kh ~ 5 x 1CT 3 s. Note that 
this is shorter than the escape time, r/v r — 0.1 s, for 
r = 100r g and v r = 0.1c. That is, there is ample time 
for the K-H instability to grow before the material flows 
outward. 

The mass-accretion rates as a function of the radius are 
displayed in Figure for the case of mi nput = 1000 and 
Z = Zq. The solid, dotted, and dashed curves indicate 
the m profiles at elapsed times of t = 10 s, 30 s, and 
50 s, respectively. Here, the mass-accretion rate at each 
radius is evaluated by 

c 2 r 

rh r = -— 2TTr 2 p(r,9)mm[0,v r (r,9)]smede. (28) 

J 

It is found that the mass-accretion rates are not constant 
in the radial direction but decrease inward, as the flow 
approaches the black hole. Roughly, we find rh r oc r. 
In addition, we find that the radial m r profile remains 
nearly the same after 7s, from which we can conclude 
that the flow is in a quasi-steady state. This rh r change 
is caused by a cooperation between a wind mass loss 
around the disk surface and the circular motion deep 
inside the disk. Note that accretion rates are assumed to 
be constant in space in the slim-disk model formulation. 

The multi-dimensional numerical simulations of RIAFs 
have shown that the accretion disks are convectively un- 
stable, thereby the circular motion being driven within 
the flow (e.g., Igumenshchev & Abramowicz 1999; Stone, 
Pringle, & Begelman 1999; McKinney & Gammie 2002). 
The convection might cooperate with the K-H instabil- 
ity in generation of the circular motion as well as the 
cavities found in our simulations. The simulations of 
RIAFs have also revealed that rh r increases with radius 
as m r oc 7- 3 / 4_1 ; which is similar to our value, and was 
attributed to the circular motion as well as the bipo- 




FlG. 6. — The mass-accretion rates as functions of the radius 
at t =10 s (by the solid line), 30 s (dotted line), and 50 s (dashed 
curve). The adopted parameters are m; nput = 1000 and Z = Zq. 



lar outflows (Stone, Pringle, & Begelman 1999, see also 
Narayan, Igumenshchev, & Abramowicz 2000 for self- 
similar solution). Here, we need to emphasize that al- 
though the flow structures look similar at first glance, 
the basic physical processes are distinct from those of 
RIAF simulations. That is, entropy is carried by fluid in 
the RIAF, whereas it is mostly by photons in the present 
case. In addition, since the RIAF simulations do not ba- 
sically take into account of the radiative cooling, they 
would overestimate the driving force of the outflows. 

4.2. Photon-trapping 

The photon-trapping characterizes the super-critical 
accretion flows. It works to reduce the energy conver- 
sion efficiency, L/Mc 2 (see, e.g. Paper I; Ohsuga et al. 
2003). As a result, the flow luminosity becomes insensi- 
tive to the mass-accretion rate when m 3> 1. However, 
simple stream lines were assumed in the previous stud- 
ies on super-critical flows, including the slim-disk model. 
Even if we do not consider possible multi-dimensional gas 
motions, photon trapping works to some degree, leading 
to a reduction in the energy conversion efficiency. Such 
multi-dimensional gas motion and associated reaction of 
the radiation can be calculated in the present RHD sim- 
ulations, since we have calculated the radiation energy 
transport, being coupled with the gas dynamics. 

Figured plots the luminosity as a function of the mass- 
accretion rate onto the black hole. The filled squares, 
circles, and triangles indicate the results for Z = WZq, 
Zq, and 0, respectively, with different mass-input rates, 
Minput = 300, 1000, and 3000 from the left to the right. 
We also indicated in the same figure the luminosity cal- 
culated based on the slim-disk model (Watarai et al. 
2000) and the one based on Model A of Paper I with 
dashed and dotted curves, respectively. In Paper I, a 
simple model for the accretion flow is employed, and the 
luminosity with carefully taking account of the photon- 
trapping is calculated by solving energy transport inside 
the accretion flows. [More precisely, the luminosity plot- 
ted in Figure is the corrected one, which fixes initial 
small errors in their Model A (see Figure 1 in Paper I).] 
Here, it should be stressed that the mass-accretion rate 
onto the black hole is not input parameter but is calcu- 
lated dynamically in the present simulations, although it 
was a parameter in both of Model A of Paper I and the 
slim-disk model. The resultant m profile was shown in 
Figure HO 

It is evident in Figure that the calculated luminos- 
ity agrees more with Model A of Paper I, rather than 
the slim-disk model, in all the cases. This proves that 
the two-dimensional effects of photon-trapping is really 
significant in the super-critical accretion flows. This re- 
sult is, in a sense, surprising. In Model A in Paper 1, 
we assumed that the viscous heating occurs only in the 
vicinity of the equatorial plane, although the gas might 
be heated up at the high altitude. Since the photons 
emitted at deep inside the disk tends to be more ef- 
fectively trapped in the flow, we anticipated that the 
photon-trapping effects would be reduced in the realistic 
situation, compared with Model A. 

We also argued in Paper I that photon trapping effects 
may be attenuated by the presence of large-scale circu- 
lation motion, which could help photon diffusion motion 
and thus considerably reduce photon traveling time to 
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Fig. 7. — The luminosities as functions of the mass-accretion 
rate onto the black hole. The symbols of triangles, circles, and 
squares indicate the metalicity of 0, Zq, and 10 Zq, respectively. 
The normalized accretion rate of each symbol is m; nput = 10, 100, 
and 1000 from the lower-left to the upper right, respectively. For 
comparison, the same but based on the slim-disk model and that 
based on Model A of Paper I are plotted by the dotted and dashed 
curves, respectively. The present results agree more with the latter 
than the former. 

the surface of the flow. However, our current results 
show significant photon trapping effects even when we 
explicitly include complex flow motions. 

Figure [7| also shows that the mass-accretion rate onto 
the black hole increases with an increase of the metalic- 
ity (for fixed mass- input rates). In our simulations the 
gas with higher metalicity has larger absorption opacities 
so that the gas energy can be more effectively converted 
into the radiation energy, yielding smaller gas pressure 
than in the metal-poor case. The gas could be effectively 
blown away by the strong gas pressure. However, there 
is a counter effect. The gas with large absorption opac- 
ity enhances the radiation energy. Enhanced radiation 
flux and large opacity can drive strong radiation-pressure 
driven outflows. The physical reason will be investigated 
in future. 

Let us see more explicitly how significant photon trap- 
ping is. Figure |H1 compares the distributions of the radial 
component of the radiative flux in the comoving frame 
(upper panel) and that in the inertial frame (lower panel) 
at t = 10 s. Other parameters are the same; the mass- 
accretion rate is TOi nput = 1000 and the metalicity is 
Z = Zq. The former radiative flux (Fq) is roughly pro- 
portional to the radial gradient of the radiation energy 
distribution. On the other hand, the latter flux (F r ) in- 
cludes the advective transport of radiation energy (v r Eo) 
in addition to the former flux; namely, we approximately 
have F r ~ F r + v r E a . In other words, differences be- 
tween two panels represent how significant photon trap- 
ping is. In fact, we see a significant difference between 
the two panels; whereas the blue area (indicating large 
radiation flux) is restricted to the vicinity of the black 
hole in the upper panel, it is more widely spread over 




Fig. 8. — The 2D distribution of the radial component of the 
radiation energy flux in the comoving frame (upper) and that in 
the inertial frame (lower), respectively. The adopted parameters 
are mj nput = 1000 and Z = Zq. The difference between them 
indicate the contribution of advective energy transport, v r Eo (see 
text). 

the entire disk region in the lower panel. This is a direct 
manifestation of the photon-trapping effects. 

4.3. Collimated Emission 

Since the super-critical accretion flows are geometri- 
cally and optically thick, the observed images and lumi- 



11 




-400 -300 



X/t c 



300 



400 




-400 -200 





-400 -200 



Fig. 9. — Surface temperature images for the model for different viewing angles: cosi = 1/8, 3/8, 5/8, and 7/8 in the upper-left, upper- 
right, lower-left, and the lower-right panels, respectively. The origin of the Cartesian coordinate on the observer's screen, (X, Y) = (0,0), 
is set at the location of the black hole. The black hole at the center The adopted parameters are rh; nput = 1000 and Z = Zq ■ 



nosity should strongly depends on the viewing (inclina- 
tion) angle. We calculate the intensity map with using 
the monochromatic radiation transfer equation, 

3 

l-VI = pT 3 (l 



v ■ l 



kB 



C p<7T 

Air m-o 



E - xh 



(29) 



where I is the specific intensity, V = (1 — v 2 /c 2 ) 1 ' 2 
is the Lorentz factor, I = (sin 9 cos $, sin sin<f>, cos 0) 



andZ = T-^l-v-l/cy^l-iT/cJv+ir-l/v^il-v^} 
are the directional cosine in the inertial and comoving 
frame, respectively, with and $ being the azimuthal 
and polar angles of radiation propagation in the inertial 
frame. Here, we assumed isotropic scattering. 

The calculated effective temperature (T c g) maps at 
t = 50 s are shown in Figure El for various viewing an- 
gles: cosi = 1/8, 3/8, 5/8, and 7/8. The other parame- 
ters are mi npu t = 1000 and Z — Zq. On the observer's 
screen, the black hole lies at the center of the Cartesian 
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coordinate (X, Y) . Since we use the spherical computa- 
tional domain with the size of 500r g in the present study, 
the blue four corners of the maps indicate outside of the 
domain. In the face-on view (see the upper left panel; 
cosi = 7/8), the effective temperature exceeds 3 x 10 6 K 
within 100r g and amounts ~ 10 7 K in the central region. 
Such a high temperature region disappears in the upper 
right panel (the case with cosi = 3/8). In this panel, 
the most luminous region is found on the upper side (not 
at the center) and is elongated in the vertical direction. 
These are caused by an occultation of the innermost part 
of the flow by the outer parts (see Fukue 2000, Watarai 
et al. 2005 for self-occultation effects based on the slim- 
disk model). Since the accretion flow is both geometri- 
cally and optically thick, the innermost region cannot be 
seen for large inclination angles, i. 

The emission from the super-critical accretion hows is 
mildly collimated due to the same reason. Figure ITU1 
represents the viewing angle-dependence of the isotropic 
luminosity (normalized by bolometric luminosity) for 
Z = Zq (circles) and Z = 10-Z© (squares), respectively. 
The other parameters are t = 50 s and rhi npu t = 1000. 
Here, the isotropic luminosity is calculated by assuming 
isotropic radiation field, L(i) — 4irD 2 I (i) £, with D be- 
ing the distance from an observer and £ being a revised 
factor. The luminosity calculated by solving radiation 
transfer equation does not always coincide with that eval- 
uated by the FLD method. Here, we assume ( = 1.7 to 
fit these luminosities. As expected, the isotropic lumi- 
nosity is quite sensitive to the viewing angle, indicating 
that the emission from the flow is mildly collimated. If 
the flow shape would be flat and (geometrically) thin, 
and if no collimation would occur, the observed luminos- 
ity should vary along the cosine curve, as is indicated by 
the dotted curve. 

We also find that the angle dependence of the isotropic 
luminosity is more enhanced in large-metalicity case of 
Z = IOZq than that of Z = Zq. This is because of dif- 
ferent density distribution. As have already mentioned, 
the mass-outflow rate is smaller in the high metalicity 
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Fig. 10. — Isotropic luminosities as functions of the viewing 
angle, i, for the cases with Z = Zq and 10 Zq. The dashed line 
represents the curve of cos i; i.e., the variation expected if the flow 
shape is flat and geometrically thin (infinitesimally). 



case than in the low metalicity case. As a result, the 
density contrast between the disk region and the out- 
flow region above the disk gets larger with an increase 
of the metalicity. Accordingly, the hot innermost region 
of the high metal flows are easier to be observed owing 
to smaller optical depth in the outflow region for small 
viewing angles. 

The super-critical accretion flows would be identified 
as very high L/Le objects in the face-on case, since the 
emission is mildly collimated in the polar direction and 
the bolometric luminosity itself exceeds the Eddington 
luminosity. The former effect is enhanced in the case of 
high metalicity (large absorption opacity). It is also true 
that super-critical flow may be identified as fairly sub- 
critical source, if the viewing angle is large (Watarai et 
al. 2005). 

5. DISCUSSION 
5.1. Comparison with the slim-disk model 
5.1.1. Flow structure 

Here, we directly compare the structure of the super- 
critical accretion flows based on our 2D simulations with 
that of the slim-disk model. Figure ^] represents the 
time-average density profile (top panel) , radial and rota- 
tion velocity profiles (middle panel), and radiative tem- 
perature profile (bottom panel) on the equatorial plane 
in the quasi-steady accretion phase of t = 40 — 50 s. The 
adopted parameters are mi npu t = 1000 and Z = Zq. The 
resultant profiles of the one-dimensional numerical study 
for the slim-disk model are represented by the dashed 
curves (Watarai et al. 2000). As have already mentioned 
in section 4.1, the injected matter falls with the free-fall 
velocity outside the Keplerian radius (shaded region) . A 
viscous accretion disk forms in the inner part (white re- 
gion). The free-falling matter penetrates slightly within 
the Keplerian radius, thus a separation radius appears 
around 70 r g . We focus on the viscous accretion regime 
in this subsection. 

The density profile is shown by the solid curve in the 
top panel. The density peak around 70r g is made by ac- 
cumulation of the injected matter, where occurs because 
of the centrifugal barrier. We compare the density pro- 
file in the viscous accretion regime with that of numerical 
and self-similar solutions of the slim-disk model. In this 
panel, it is found that the slope of the density profile 
of our simulations is steeper than that of the numeri- 
cal solution of the slim-disk model. We also find that 
the slope is close to that of the self-similar solution of 
the slim disk, p oc R-^l 2 (Watarai et al. 1999, see also 
Spruit et al. 1987), whereas the profile becomes flatter 
in the vicinity of the black hole (R < 5r g ). Here, we note 
that the slope in our simulations is steeper than that of 
the RIAFs. The RIAF simulations showed the density 
profile to be p oc r -1 / 2 (McKinney & Gammie 2002). 

In the middle panel, the radial and rotation velocities 
are plotted by the filled circles and the solid curve. The 
dotted curve indicates the Keplerian velocity, wkep = 
V GMR/(R — r s ). It is found that the rotation velocity 
nicely agrees with the Keplerian velocity as well as the 
result of the numerical solution of the slim disk. How- 
ever, the radial velocity profile is not smooth and largely 
deviates from that of the slim-disk value, whose slope 
agrees with that of the Keplerian velocity (dotted curve) 
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Fig. 11. — The density profile (top), radial and rotation velocity 
profiles (middle), and radiation temperature profile (bottom) at 
the equatorial plane, z = 0, averaged in the elapsed time between 
t = 40 s and 50 s. The adopted parameters are mj nput = 1000 
and Z = Zq. Shaded region is the free-falling region, where the 
injected matter falls at the free-fall velocity, whereas the viscous 
accretion disk forms around the black hole, r < 70r g . The dashed 
curves indicate the profiles of one-dimensional numerical study of 
the slim-disk model (Watarai et al. 2000). The thin solid lines 
in the top and bottom panel means the slope of the self-similar 
solution of the slim-disk model. The dotted curve in the middle 
panel is the Keplerian velocity. 

in the free-fall regime. The radial velocity remarkably 
varies at R = 5 — 70 r g , and the radial velocity some- 
times becomes negative, around R — 6, 20, and 60 r g . 
Such a complex v r profile is formed by the prominent 
circular motions around the equatorial plane (see Figure 
0). 

The bottom panel shows the radiation temperature 
profile, where the radiation temperature is defined as 
T r = (E/a) 1 / 4 with a being the radiation constant. As 
shown in this panel, the profile nicely agrees with that 
of the one-dimensional numerical solution and is flat- 
ter than that of the self-similar solution of the slim-disk 
model, T r oc i?~ 5 / 8 . 

To sum up, we can thus conclude that the quasi-steady 
structure of the super-critical accretion flows are appar- 
ently similar to but, precisely speaking, deviates from 
the prediction of the slim-disk model. 

5.1.2. Effective temperature 

Next, we represent again the effective temperature pro- 
file in order to directly compare with the prediction of the 
slim-disk model, although the 2D images of the effective 
temperature have already shown in Figure Figure IT2"1 
represents the effective temperatures at Y = as a func- 
tion of X for i = and tt/6, where X and Y are the hor- 
izontal and vertical coordinates on the observer's screen 
(see Figure EJ- The other parameters are m input = 1000, 
Z = Zq, and t — 50 s. It is found that the slope of the 
effective temperature profile for i — becomes flatter 
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Fig. 12. — The effective temperature profiles in the X-directions 
at Y = for the case with i = (solid curve) and 7r/6 (dotted 
curve). Here, (X,Y) is the Cartesian coordinate on the observer's 
screen, and the black hole lies at the center. The other parameters 
are m; nput = 1000, Z = Zq, and t = 50 s. The dashed curve 
indicates the result of the one-dimensional numerical study of the 
slim-disk model (Watarai et al. 2000), and the thin solid line means 
the slope of the self-similar solution of the slim-disk model. 

than that of the slim-disk model, T c ff cx i?" 1 / 2 , within 
30r g . In contrast, we found that the profile is consistent 
with the slim-disk model in the outer region. As have 
already mentioned, the effective temperature as well as 
the luminosity is quite sensitive to the viewing angle, 
since the accretion flow is both geometrically and op- 
tically thick. We found that the effective temperature 
profile for z = tt/6 is almost flat. 

The slim-disk model succeeds in reproducing the ob- 
served behavior of Ultra-luminous X-ray sources (ULXs) 
and narrow- line Seyfert 1 galaxies (NLSls), which are 
thought to be candidates of near- or super-critical flows 
(Watarai, Mizuno, & Mineshige 2001; Mineshige et al. 
2000; Kawaguchi 2003). The SEDs are produced based 
on the blackbody emission with the effective temperature 
coupled with the some modification. We, here, note that 
the temperature profile obtained by our current simu- 
lations show some deviations from that of the slim-disk 
model. 

5.2. Blueshifted absorption by outflow matter 

One characteristic feature of the super-critical accre- 
tion flow is its generation of the radiation-pressure driven 
wind. Our simulations show that a large quantity of gas 
is blown away by the strong radiation pressure above 
and below the accretion disk. Such outflow material, 
which is thought to be highly ionized by the radiation 
from the accretion flow, would absorb the continuum 
emission. Therefore, the super-critical accretion objects 
have the blueshifted absorption lines by the highly ion- 
ized ions. Such blueshifted absorption in high-ionization 
lines was observed in the UV spectra of broad absorption 
line quasars (Weymann, Carswell, & Smith 1981; Becker 
et al. 2000). Moreover, Pounds et al. (2003a, 2003b) 
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and Reeves, O'Brien, & Ward (2003) reported by the 
blueshifted X-ray absorption lines that the highly ionized 
matter with outflow velocity of order ~ 0.1c in quasars, 
and the flow is likely to be optically thick to the electron 
scattering within r ~ 100r g . In our simulations, the typ- 
ical outflow velocity is 0.1c and the optical depth of the 
outflow regions is around 1.2(p/10 _8 gcm _3 )(r/100r g ) 
(This estimated value does not depend on the black- 
hole mass because of p cx M and r oc M _1 .) These 
are remarkably close to the observed results, and there- 
fore, our simulations basically account for the origins of 
high-velocity outflows. Detailed comparison with the ob- 
servations is left as future work. 

5.3. Growth of supermassive black holes 

We reveal by 2D-RHD simulations that the black hole 
can swallow the gas at the super-critical accretion rate, 
although the luminosity exceeds the Eddington luminos- 
ity. This result implies that the black hole can rapidly 
grow, and the growth timescale is given by M/M = 
4.5 x 10 6 (m/100) _1 yr. Thus, if the supermassive black 
holes in the galactic nuclei are built up by the super- 
critical accretion, such a growing phase would be quite 
short because of a short growth timescale. 

Umemura (2001) suggested the radiation- 
hydrodynamical model for the formation of the 
supermassive black holes, in which the supermassive 
black holes grow in the ultraluminous infrared galaxies 
(ULIRGs) due to the mass accretion caused by the radi- 
ation drag. Based on this model, Kawakatu, Umemura, 
& Mori (2003) also claim that proto-quasars, which 
have relatively small black hole and show only narrow 
emission-lines, appear just after the ULIRG phase, 
t ~ 10 8 ~ 9 yr. However, such proto-quasar phase might 
not be observed if the black hole grows at super-critical 
accretion rate in ULIRG phase. A seed black hole 
in the ULIRG can evolve to a supermassive black 
hole via the super-critical accretion by the end of the 
ULIRG phase, since the mass-accretion rate from the 
circumnuclear regions can exceed the critical rate due 
to the effective radiation drag. Thus, we presumably 
observe the grown-up quasars after the ULIRG phase, 
and proto-quasars are obscured by huge amount of the 
dust in the ULIRG. 

As shown in Figure EI our simulations show that the 
black hole can effectively grow in the case that the ac- 
creting gas is metal rich. This implies that super-critical 
accretion flows tend to emerge in the metal rich objects 
like star burst galaxies or star- forming regions. This ten- 
dency is also consistent with the observed results that 
the NLSls, which is though to be near- or super-critical 
accretion objects, are metal rich objects (Nagao et al. 
2002; Shemmer & Netzer 2002; Shemmer et al. 2004). 
It is proposed that NLSls accrete at super-critical rates 
based on the study of the optical band (Kawaguchi 2003; 
Collin & Kawaguchi 2004). 

The growth of the black hole may gradually slow down, 
since the density as well as the absorption opacity of 
the super-critical accretion flows are small around the 
massive black hole, if the normalized mass-accretion rate 
does not change so much. Observational constraints are 
proposed by Yu & Tremaine (2002), by which most of 
luminous quasars must be sub-critical accretion phase. 
They showed that the luminous quasars have the energy 



conversion efficiency of L/Mc 2 > 0.1, based on the study 
of the local black hole density and the luminosity func- 
tion of the quasars. Since this efficiency is much smaller 
than 0.1 in the super-critical accretion flows, most of lu- 
minous quasars must be sub-critical accretion phase (see 
also Soltan 1982). 

5.4. Future work 
5.4.1. Spectral energy distribution 

Throughout the present study, we use the frequency- 
integrated energy equation of the radiation [see equation 
(jHJ)]. By solving the monochromatic radiation transfer 
equation, we can calculate the effective temperature pro- 
files and show them in Figure However, the emergent 
SEDs of the super-critical accretion flows are not a simple 
superposition of blackbody spectra with various effective 
temperatures (Ohsuga et al. 2003). The photons gener- 
ated deep inside the disk are difficult to reach the disk 
surface and thus moves downward with gas. Further- 
more, most of photons generated in the vicinity of the 
black hole will be immediately swallowed by the black 
hole and cannot contribute to the emergent SED. Fre- 
quency dependent radiation-hydrodynamical simulations 
are necessary to resolve this issue. 

According to the current simulations the hot outflow 
with temperature of 10 9 ~ 10 K appears above the disk. 
The Compton y-parameter of this outflow region is com- 
parable to or larger than unity, meaning that the inverse 
Compton scattering cannot be ignored. Seed soft pho- 
tons emitted from the disk region should be Compton 
up-scattered, producing high-energy, non-thermal emis- 
sion which contributes to the emergent SED. The corona 
above the super-critical accretion disk has been claimed 
to fit the observed SEDs of NLSls and very-high state of 
BHCs (Wang & Netzer 2003; Kubota & Done 2004), al- 
though the formation mechanism of the corona is poorly 
known. The hot outflow shown in our simulations might 
resolve this issue. 

Here, we should note that the Comptonization process 
is not included in our simulations. If the Compton cool- 
ing were effective, the hot outflow might be cooled to 
some extent. The bulk Comptonization may also con- 
tributes the production of the high-energy photons be- 
cause of large radial velocities, v r > 0.1c, near the black 
hole. The detailed study of the SED with Comptoniza- 
tion is, however, beyond the scope of this paper and 
should be explored in future work. 

5.4.2. Viscosity 

We assumed that only the rc^-componcnt of the viscous 
stress tensor is non-zero while all the other components 
vanish, because the r(p-component plays an essential role 
of in the angular momentum transport within the disk. 
If the r9- as well as (^-component were non zero, the 
structure near the disk surface might change. Figure |31 
clearly shows abrupt velocity changes across the bound- 
ary between the disk and the outflow regions, which is 
thought to promote the K-H instability. The growth 
of this instability might be suppressed by the rO- and 
^-components of viscosity. These components might 
also partially suppress the circular motion in the disk. 
Stone et al (1999) discovered by the simulations of RI- 
AFs that the flow structure differs from that given by 
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Igumenshchev & Abramowicz (1999). It implies that the 
inclusion of the components of the viscous stress tensor 
suppresses the convection, since the r#-component is ex- 
cluded in Stone et al. (1999). This will be examined in 
future study. 

More importantly, we need coupled MHD and RHD 
simulations, since the source of disk viscosity is likely to 
be of magnetic origins (Hawley, Balbus, & Stone 2001; 
Machida, Matsumoto, & Mineshige 2001; Balbus 2003 
for a review) . Recently, local radiation MHD simulations 
of the accretion flow have been performed by Turner et 
al. (2003) and Turner (2004). We definitely need global 
simulations to include global field effects. 

5.4.3. The FLD approximation 

Finally, we need to comment on the FLD approxima- 
tion. It is known that the FLD approximation does not 
always give good results for the regions with moderate 
optical thickness. The FLD approximation is thought 
to be valid in the disk region, which has a large optical 
depth, but the spherical shape of the radiation energy 
density distribution above the disk might be artificial. 

In addition, the radiation drag force cannot be treated 
in the FLD approximation, whereas it might play an im- 
portant role for the transport of the angular momentum 
within the outflow region. Since the timescale of the an- 
gular momentum transport via the radiation drag, which 
is given by 




(with M g being the mass within the outflow region), is 
about ten times longer than the escape time, r/v r — 
0.2s(r/200r g )(u r /0.1c)" 1 (M/10M Q ), roughly ten per- 
cent of the angular momentum could be extracted in the 
outflow region. [The exact expressions for the radiation 
drag are found in the literature (Umemura, Fukue, & Mi- 
neshige 1997; Fukue, Umemura, Mineshige 1997; Ohsuga 
et al. 1999).] 

Radiation drag arises where there exists a large veloc- 
ity difference between radiation sources and the irradi- 
ated matter. In the FLD approximation, however, the 
radiation flux is determined solely by the local gradient 
of the radiation energy, and photon trajectories cannot 
be properly considered. Thus, the FLD approximation 
can not treat the radiation drag force, in principle. It 
would be better to solve the radiation transfer equations 
without using this approximation. 

Bcgelman (2002) suggested that strong density inho- 
mogeneities on scales much smaller than the disk scale 
height is formed due to the photon bubble instability 
in the radiation pressure-dominated accretion disks, and 
the disk can remain geometrically thin even as the the 
maximum luminosity exceeds Eddington luminosity by 
a factor of one hundred. However, the FLD approxima- 
tion is not suitable method for investigating the radia- 
tion fields in such inhomogeneous structure. The detailed 
study of the photon bubble instability should be explored 
in future work. 



6. CONCLUSIONS 

By performing the 2D-RHD simulations, we, for the 
first time, investigate the quasi-steady structure of the 
super-critical accretion flows around the black hole with 
particular attention being paid on the photon-trapping 
effects. We have obtained several new findings: 

(1) The quasi-steady structure of the super-critical flow 
is divided into two parts: the disk region (with mostly 
inflow) and the outflow region above the disk. The two 
regions are separated by a sharp density jump. The gas 
outflow driven by the strong radiation pressure is pro- 
duced around the rotation axis. Further, there exists 
velocity shears at the boundary, which causes K-H in- 
stability around the disk surface, producing the patchy 
structure as well as the circular motion within the disk 
region. Convection may also be responsible for such in- 
homogeneous structure. 

(2) The photon-trapping plays an important role in 
the super-critical accretion regime. The advective energy 
transport is substantial, and the large amount of photons 
generated inside the disk is swallowed by the black hole 
without radiated away. 

(3) Our 2D-RHD model shows some differences from 
the slim-disk model. The slim-disk model assumes a 
simple convergence flow, while our simulations revealed 
rather complex gas motion and structure. We also found 
that the mass-accretion rate is not constant in space but 
decreases as matter accretes, roughly as m oc r, as a 
result of wind mass loss and circular motion. The calcu- 
lated luminosity of the flows agrees more with the predic- 
tion of Paper I rather than that of the slim-disk model. 

(4) The emission of the super-critical accretion flows is 
moderately collimated. The apparent luminosity could 
become more than ten times larger than the Edding- 
ton luminosity. The super-critical accretion flows would 
be identified as high L/Le objects in the face-on view, 
but not, if the viewing angle is large, for which self- 
occultation tends to reduce the total luminosity and the 
maximum flow temperature. 

(5) The mass-accretion rate increases with increase of 
the absorption opacity (metalicity) of the accreting mat- 
ter. It implies that the black hole tends to grow up faster 
in the metal rich regions as in starburst galaxies or star- 
forming regions. In addition, the growth of the black 
hole may gradually slow down, since the density as well 
as the absorption opacity of the super-critical accretion 
flows are small around the massive black hole, if the nor- 
malized mass- accretion rate is kept constant. 
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